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On  Calculating  Analytic  Centers 

A.  A.  Goldstein* 


This  note  was  motivated  by  papers  of  Renegar  and  Shub(88)  and  by  Ye(89).  We  apply 
Smale's(86)  estimates  at  one  point  for  Newton's  method  to  the  problem  of  finding  the 
analytic  center  of  a  polytope.  The  method  converges  globally  in  the  appropriate  norm. 
The  ideas  are  then  applied  to  obtain  a  possible  benchmark  for  path  following  methods. 

When  Smale's  method  is  tractable  its  power  stems  not  only  from  the  fact  that  the  in- 
formation is  concentrated  at  one  point.  There  are  2  norms  to  estimate,  not  3  as  in  the 
Kantorovich  estimate.  Moreover  no  estimate  of  the  inverse  of  the  derivative  operator  by 
itself  is  needed.  The  need  for  the  norm  of  the  inverse  by  itself  often  makes  for  coarse 
estimates. 

1.   Setting 

Let  A  denote  an  m  by  n  orthonormal  matrix  of  rank  n  and  b  an  m  by  1  matrix. 
We  assume  that  m  >  n.  Denote  by  e  a  m  by  1  matrix  whose  components  are  all  ones. 
Transposes  of  matrices  will  be  denoted  by  an  asterisk,  rows  of  a  matrix  by  superscripts, 
and  columns  by  subscripts.  The  euclidean  space  of  real  m-tuples  will  be  denoted  by  Em. 
If  u  E  Em  we  mean  by  diag(u)  the  diagonal  matrix  with  entries  utt.  The  dot  product 
corresponding  to  the  usual  norm  will  be  denoted  by  ,  ].  The  usual  norm  will  be  written 
as  ||  ||2.  En  will  be  also  be  renormed  under  a  dot  product  that  will  be  denoted  by  <  ,  >. 
The  norm  arising  from  this  dot  product  will  be  written  as  ||  ||.  Let  P  be  a  polytope  with 
non-empty  interior  given  by  the  inequalities 

b  -  Ax    >    0. 

Given  Xq  in  the  interior  of  P  and  e  >  0,  we  seek  the  analytic  center  £  of  P  to  within  a 
tolerance  of  e.  Let  Rt(x)  =  bl  —  Alx  . 

Claim  1.    Let  N  be  the  smallest  integer  exceeding 


log2(4.95m4  maxRt(x0))    +   log2  (  - 


1  +  log2 

Then  if  N  steps  of  the  Newton  sequence  are  generated  using  the  gradient  of  the  potential 
function  below 
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||* JV  -  CII2    <    e 
The  proof  of  the  claim  depends  on  the  following  ingredients. 

2.  Some  ingredients. 

The  potential  of  P  is  defined  by  the  expression 


?r(P)  =mai]  [Rt(x)  :x  £  P} 


1  =  ] 


The  maximum  is  achieved  at  an  unique  point  called  the  analytic  center  of  P.     (Ye  87  ). 
We  shall  find  this  point  by  seeking  a  zero  for  the  gradient  of  the  logarithmic  potential 


4>(x)  =  Y/^g(Rl(x) 


i=i 


Let  D(x)  =  diag(l/Ri(x)),  thus  Dtl{x)  =  1/Ri(x). 

We  apply  Newton's  method  to  the  gradient  of  (f>  which  we  denote  by  F.  F(x)  may  be 
represented  by  the  matrix  A*D(x)e,  and  F(x)  G  En.  The  kth  Frechet  differential  of  F  at 
x  can  be  identified  with  a  multi-linear  mapping  from  (E^j)k  to  En.  A  representation  of 
these  differentials  as  matrices  follows. 

F'(a-)/ii  =  -A*D2(x)diag(Ah1)e  =  -A*D(x)D(x)Ah1 

F"(x,huh2)  =  -2\A*  D{x)D2{x)  diag(Ah2)  Ahx 

F'"{x,huh2M)  =  -3\A*D{x)D3(x)diag(Ah3)diag{Ah2)Ah1 

and 
F(k)(xyhuh2,...,hk)  =  -k\A*  D{x)Dk{x)  diag{Ahk), ...,  diag{Ah2)Ah1 

=  -k\A*D(x)Q(x,hu...,hk) 
Here 

||Q(.r)||2  =  sup{\\Q(x,hu-,hk)h  :  INI2  =  IMk  =,  -,=  \\hkh  =  1}    <   1 


Theorem  1.  (SmaJe  86  )  Assume  F  is  an  analytic  map  between  real  Banach  spaces  X  and 

Y,  that  is  the  Frechet  derivatives  F^k\x)  exist  for  all  x  €  X  and  k=l,2,3, Given  x0  G 

X,  assume  that  the  inverse  of  F'(x)  which  we  denote  by  F'_1(x)  exists.  Set 

f3(x0)  =  \\F'_1(x0)F(x0)\\         and 

7(*o)  =  supi[\\^F'_1(xo)F^(xo)\\^  :  k  >  2 J 

If 

P(xo)f(x0)    <    .130707 

then  xo  is  an  approximate  root  of  F.  That  is,  the  Newton  sequence 

Xk+i  =  Xk  -  F'_1(xk)F(xk) 
is  well  defined  and  {x^}  converges  to  say  £,  a  root  of  F  at  the  rate: 

Moreover 

11**- £11    <    7-(\f'lf3(x0)  (A) 

3.  Proof  of  Claim  1. 

Assume  x0  is  given  in  P(M).  The  matrix 

P(x0)  =  D(x0)A(A*D(x0)D(x0)A)-1A*D(x0) 

maps  each  point  in  Em  to  its  closest  point  in  the  range  of  the  matrix  D{xq)A  .    Hence 
||P(a-o)||2  =  1.  Werenorm  En  by 

11*11  =  \\DcAx\\2 

Here  Dq  —  CD(xG)  with  C  =  1/S^m^ .  With  this  definition  we  get: 

P(x0)  =  C\\P(x0)e\\2  =  m*/si 

Also 

1(xo)<Csup(\\P(x0)\\^T)sup(\\QkAh1\\^)    <    C 

Thus         [3{x0)~f{x0)    <    -    <    .130707 
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Hence  by  Smale's  theorem  the  sequence  generated  by  the  Newton  algorithm  converges  to 
the  analytical  center  £  with  a  rate  given  by  (A)  in  Theorem  3.1  above. 

Since  <  x,  x  >—  [DcAx,  DqAx]    >    C2  \\x\\\ / max Rx(x 0)2 ,  then 

||x||2    <    CRi{xQ)\\x\\ 

Now  choose  N  so  that 

CRi(xQ)\\xN-Z\\    <    e 

4.  Application  to  programming 

By  a  theorem  of  Ye  (89),  if  one  of  the  hyperplanes  of  P  is  translated  to  pass  thru  £  then 
the  resulting  polytope  P+  satisfies 

tt(P+)  1 

— <    - 

tt(P)      -    e 


Consider  the  following  algorithm  for  linear  inequalities.  We  wish  to  solve  the  system 
b  —  Ax  >  0  if  this  is  possible.  Given  an  arbitrary  xq  choose  M  so  that  b  -f  M  —  Ax  >  0. 
Find  the  center  of  this  polytope  P(M).  Take  the  smallest  component  of  .R(£),  say  Rqd)- 
Begin  anew  with  the  polytope  P(M  —  Rq(£)).  This  algorithm  has  a  worst  case  iteration 
count  of  0(777.)  times  our  cost  of  getting  to  the  center. 

For  linear  programming  let  the  polytope  P  be  given  by  b  -  Ax  >  0  and  P(M)  the  polytope 
define  by  the  inequalities  for  P  together  with  the  inequality  M  -  [c  ,  x]  >  0.  We  seek  the 
smallest  M  for  which  P(M)  is  non-empty.  We  first  find  the  center  £  of  the  polytope  P. 
We  then  find  the  intersection  of  the  ray  {x  =  £  —  tc  :  t  >  0}  with  P  Translate  the  cost 
hyperplane  to  pass  thru  this  point.  Then  find  the  center  of  the  new  polytope  P(M). 

5.  Benchmark 

We  now  consider  the  possibility  of  starting  from  a  point  in  a  polytope  P(M)  and  moving 
to  the  center  of  a  neighboring  polytope  P(M  —  l/2y/m)  by  Newton  steps. 


Assume  that  at  (x0 ,  M0),  Ri(x,  M)  —  bt  +  M0  —  A%x    >   0.  We  seek  a  point  (x1:  Mi)  such 
that 

=    0,  1  <j  <  n  (la) 


a<A(x,M) 


dx} 
d<Kx,M)         ^(n.MQ  .         MM 

~M m—  =  °      {lb) 
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and  such  that 


Riibi  +  Mi  -  Axxx)    >   0         (2) 

Let  Mi  =  M0  —  l/2y/m.  Assume  that  the  value  of  X\  is  well  defined  and  given.  Otherwise 
P(Mi)  is  empty  and  Mq  is  within  l/2y/m  of  M*  the  optimal  value  of  M.  We  show  that 
(aro,Mo)  is  an  "approximate  root"  for  system  (I). 

Remark  The  matrix  (A  e)  has  rank  n+1. 

Proof  Because  of  our  boundedness  assumption  on  the  polytopes,  the  system  of  inequalities 
Ax  >  0  is  inconsistent.  If  u  is  in  the  null  space  of  (A  e)  then  Au  =  -un+ie  ^  0,  a 
contradiction. 

In  matrix  notation  the  system  (1)  (after  scaling  the  second  entry)  is 


-A*D{x,M)e         -^=e*(D(x,M)e-D(xuM1)e) 


Thus  we  see  that  -  F'(x,  M)  may  be  generated  from  the  matrix 


—  g 
B  =  (A1,A2,...,An,An+1)         where  An+1  = 


2y/m 

Assume  that  (A\,  A2, ...,  An)  is  rescaled  if  necessary  so  that  ||J5||  <  1.  By  the  Remark  we 
see  that  B  has  rank  n+1.  Thus  Claim  1  holds  for  this  case  as  well.  If  we  are  satisfied  with 
a  reduction  of  1/3-^/m  this  will  happen  in  N  steps  by  the  claim  with  e  set  to  l/6>/m.  We 
have  then  the  following  result:  (not  an  algorithm  but  a  benchmark) 

Claim  2.  We  are  given  a  point  (xjt,Mjfc),.  Let  Mk+i  —  M*  -  l/2v/m.  If  P(Mk+i)  is 
not  empty,  take  £jt+i  for  its  center.  Let  the  system  (I)  be  run  with  Newtons'  method. 
Otherwise,  stop.  In  N  steps  M*  will  be  reduced  by  at  least  l/3>/m.  This  value  updates 
Mk+i  and  the  corresponding  iterate  for  x  updates  xjt+i-  Assume  the  optimal  M  say  M* 
known.  Then  the  global  Newton  process  can  be  terminated  in  no  more  than  Q  steps,  where 


Q>  3^i(M0  -M*)(l  +  log2 


log2 (4.95m4  maii?,(xo))    +   log2 (  6\/m 


At  termination  M/v  is  within  l/2y/m  of  M*  and  x ^  is  an  approximate  root  for  system   (I) 
with  M*  replacing  M\  and  £  replacing  x\,  respectively. 

A  similar  result  holds  for  linear  programming. 
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